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ABSTRACT 

Using Fourier analysis, we prove that the limiting distribution of the standardized random 
number of comparisons used by Quicksort to sort an array of n numbers has an everywhere pos- 
itive and infinitely differentiable density /, and that each derivative enjoys superpolynomial 
decay at ±00. In particular, each is bounded. Our method is sufficiently computational to 
prove, for example, that / is bounded by 16. 

AMS 2000 subject classifications. Primary 68W40; secondary 68P10, 60E05, 60E10. 

Key words and phrases. Quicksort, density, characteristic function, sorting algorithm, Fourier 
analysis, rapidly decreasing C°° function, tempered distribution, integral equation. 

Date. March 31, 2000. 



1 Research supported by NSF grant DMS-9803780, and by the Acheson J. Duncan Fund for the 
Advancement of Research in Statistics. 



1 



1 Introduction and summary 

The Quicksort algorithm of Hoare |?j is "one of the fastest, the best-known, the most 
generalized, the most completely analyzed, and the most widely used algorithms for 
sorting an array of numbers" 0. Quicksort is the standard sorting procedure in Unix 
systems, and Philippe Flajolet, a leader in the field of analysis of algorithms, has noted 
that it is among "some of the most basic algorithms — the ones that do deserve deep 
investigation" j|. Our goal in this introductory section is to review briefly some of what 
is known about the analysis of Quicksort and to summarize how this paper advances 
that analysis. 

The Quicksort algorithm for sorting an array of n numbers is extremely simple to 
describe. If n = or n = 1, there is nothing to do. If n > 2, pick a number uniformly 
at random from the given array. Compare the other numbers to it to partition the 
remaining numbers into two subarrays. Then recursively invoke Quicksort on each of 
the two subarrays. 

Let X n denote the (random) number of comparisons required (so that Xq = 0). 
Then X n satisfies the distributional recurrence relation 

X n k X Un - X + X*_ Un + n - 1, n > 1, 

where = denotes equality in law (i.e., in distribution), and where, on the right, U n is 
distributed uniformly on the set {1, . . . , n}, X*- = Xj, and 

U n ) Xq,... , X n -i; Xq,... ,X n _ 1 

are all independent. 

As is well known and quite easily established, for n > we have 

ji n := EX„ = 2(n + l)H n — An ~ 2nlnn, 

where H n := 5Zfc=i ^ s the nth harmonic number and ~ denotes asymptotic equiva- 
lence. It is also routine to compute explicitly the standard deviation of X n (see Exercise 

6.2.2-8 in ||), which turns out to be ~ n \j^ ~ l 71 " 2 - 
Consider the standardized variate 

Y n := (X n - fJL n )/n, n > 1. 

Regnier [|ll]] showed using martingale arguments that Y n — > Y in distribution, with Y 
satisfying the distributional identity 

(1.1) Y = UY+(l-U)Z + g(U)=:h Y:Z (U), 

where 



(1.2) 



g(u) := 2nlnn + 2(l - u) ln(l - u) + 1, 
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and where, on the right of = in (|1 . 1[) , U, Y, and Z are independent, with Z = Y and 



U ~ unif(0, 1). Rosier |12| showed that ( |1 . 1[ ) characterizes the limiting law C(Y), in the 
precise sense that F := C(Y) is the unique fixed point of the operator 

G = C{V) ' ^ SG := £(C/F + (1 - C/)y* + #(£/)) 

(in what should now be obvious notation) subject to 

BV = 0, VarF < oo. 

Thus it is clear that fundamental (asymptotic) probabilistic understanding of Quick- 
sort's behavior relies on fundamental understanding of the limiting distribution F. In 
this regard, Rosier [12| showed that 

(1.3) the moment generating function (mgf) of Y is everywhere finite, 

and Hennequin Q || and Rosier showed how all the moments of Y can be pumped out 
one at a time, though there is no known expression for the mgf nor for the general pth 
moment in terms of p. Tan and Hadjicostas |l5| proved that F has a density / which is 
almost everywhere positive, but their proof does not even show whether / is continuous. 

The main goal of this paper is to prove that F has a density / which is infinitely 
differentiable, and that each derivative f^ k \y) decays as y — > ±oo more rapidly than 
any power of this is our main Theorem |3.1| In particular, it follows that each 



is bounded (cf. Theorem 3.3). 



Our main tool will be Fourier analysis. We begin in Section || by showing (see 



Theorem 2.9) that the characteristic function <f) for F has rapidly decaying derivatives 
of every order. Standard arguments reviewed briefly at the outset of Section || then 
immediately carry this result over from to /. Finally, in Section 4 we will use the 
boundedness and continuity of / to establish an integral equation for / (Theorem 4.1). 
As a corollary, / is everywhere positive (Corollary 



Remark 1.1. (a) Our method is sufficiently computational that we will prove, for 
example, that / is bounded by 16. This is not sharp numerically, as Figure 4 of |15| 
strongly suggests that the maximum value of / is about 2/3. However, in future work we 
will rigorously justify (and discuss how to obtain bounds on the error in) the numerical 
computations used to obtain that figure, and the rather crude bounds on / and its 
derivatives obtained in the present paper are needed as a starting point for that more 
refined work. 

(b) Very little is known rigorously about /. For example, the figure discussed in (a) 
indicates that / is unimodal. Can this be proved? Is / in fact strongly unimodal (i.e., 
log-concave)? What can one say about changes of signs for the derivatives of /? 

(c) Knessl and Szpankowski (8| purport to prove very sharp estimates of the rates 
of decay of f(y) as y — > — oo and as y — > oo. Roughly put, they assert that the left 
tail of / decays doubly exponentially (like the tail of an extreme-value density) and 
that the right tail decays exponentially. But their "results" rely on several unproven 
assumptions. Among these, for example, is their assumption (59) that 

Ee~ Ay ~ exp(aA In A + [5\ + 7 In A + 5) as A — > 00 



3 



for some constants a(> 0), j3, 7, S. [Having assumed this, they derive the values of a, 7, 
and 6 exactly, and the value of (5 numerically.] 

2 Bounds on the limiting Quicksort characteristic function 

We will in this section prove the following result on super polynomial decay of the char- 
acteristic function of the limit variable Y . 

Theorem 2.1. For every real p > there is a smallest constant < c p < 00 such that 
the characteristic function <p(t) := Ee* <y satisfies 

(2.1) \(j)(t)\ < c p \t\- p for all t G R. 

These best possible constants c p satisfy cq = 1, C]/ 2 < 2, c 3/ / 4 < V8tt, c\ < 4ir, c 3 / 2 < 
187, c 5/2 < 103215, c 7/2 < 197102280, and the relations 

(2-2) c^<c p /^, 0< Pl < P2] 

(2.3) c p+1 < 2P +1 cl + ^Mp/(p - 1), p > 1; 

(2.4) c p < 2 p2+6p , p > 0. 

[The numerical bounds are not sharp (except in the trivial case of Co); they are the 
best that we can get without too much work, but we expect that substantial improve- 
ments are possible.] 



Proof. The basic approach is to use the fundamental relation (1.1). We will first show 



using a method of van der Corput |lj, |10| , that the characteristic function of h y:Z (U) is 
bounded by 2|£| -1//2 for each y, z. Mixing, this yields Theorem 2.1 for p = 1/2. Then 
we will use another consequence of (^J), namely, the functional equation 

(2.5) 4>(t) = [ <j>{ut) 0((1 - u)t) e it9( - u) du, t G R, 

Ju=0 

or rather its consequence 

(2-6) m)\ < T 1^)1 W(l-«)t)l du, 

Ju=0 

and obtain successive improvements in the exponent p. 

We give the details as a series of lemmas, beginning with a standard calculus estimate 



10(1 . Note that it suffices to consider t > in the proofs because (j>(—t) = 4>(t) and thus 
4>{— 1)\ = |(/>(i)|. Note also that the best constants satisfy c p = sup t>0 t p \4>{t)\ (although 



we do not know in advance of proving Theorem |2.l| that these are finite), and thus 
cj p = sup <>0 t\cj){t)\ l l p , which clearly satisfies ( |2.2| ) because \4>(t)\ < 1. 

Lemma 2.2. Suppose that a function h is twice continuously differentiable on an open 
interval (a, b) with 

h'(x) > c> and h"(x) > for x G (a, b). 
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Then 



<i ith ^ dx 



< — for all t > 0. 

~ ct 



Proof. By considering subintervals (a + e, b — e) and letting e — > 0, we may without loss 
of generality assume that h is defined and twice differentiable at the endpoints, too. 
Then, using integration by parts, we calculate 



I e ith ^ dx = - f 

Jx=a J x=a 



d c ith(x) 
dx 



1 

it 



ith(x) 



h'(x) 



dx 
h'(x) 
b 



Jth(x) j 



h'(x) 



So 



3 Uh{x) dx 



< 



1 [ / 1 



+ 



t \ \h'(b) h'(a) 

1 f f 1 1 

+ 



t \\h'{b) h'(a) 

1 fY 1 1 

+ 



+ 



+ 



+ 



t \ \h'{b) h'{a) ) \h'(a) h'(b) 



x=a 
b 

x=a 
1 



h'(x) 



dx 



1 



h'{x) 
1 



d.x 



2 2 
< — . 



□ 



th'(a) ~ ct' 

Lemma 2.3. For any real numbers y and z, the random variable hy )Z {U) defined by 
(1.1) satisfies 

| Ee ith»,.(lO| < 2\t\- 1/2 . 



Proof. We will apply Lemma 2.2, taking h to be h y>z . Observe that 



^ („) = 2 (! + —!- 

y ' \u 1 — it 



it / it(l — it) 



> 8 for n e (0, 1) 



and that 



h'v z( u ) = if and only if u = a 



1 + exp (§(y - z)) 



6(0,1). 



Let t > and 7 > 0. If in Lemma 22 we take a := a yjZ + jt l / 2 and 6 := 1, and assume 
that a < b, then note 

h'{u) = h' y z (u) = / hy jZ {x) dx > 8{u - a V)Z ) > 8jt~ 1/2 for all u G (a, 6). 



So, by Lemma 2.2, 



Jth y , z (u) du 



=Q J/ , Z +7*" 1/2 



< 2 [g r l/2l-l = -l/2_ 
~ t L ' J 4 7 



Trivially, 



so we can conclude 



< it' 



-1/2 



3 «h»,«(«) d„ 



< [(4 7 )- 1 + 7 ]f 



-1/2 



This result is trivially also true when a = a y ^ z + 7 t x / 2 > b = 1, so it holds for all 
t, 7 > 0. The optimal choice of 7 here is 7 op t = 1/2, which yields 



Jthy, Z (u) ^ 



< r 1 ' 2 for all t > 0. 



Similarly, for example by considering u 1— > h(l — u) 



< r 1 / 2 for all t > 0, 



and we conclude that the lemma holds for all t > 0, and thus for all real t. 
Lemma 2.4. For any real t, \(f>(t)\ < 2|i|~ 1 / 2 . 



□ 



Proof. Lemma 2J5 shows that 

E(t'" ly - zil ' } 

and thus 



Y,Z 



< 2\t\ 



-1/2 



Y,Z 



□ 



The preceding lemma is the case p = 1/2 of Theorem 2A. We now improve the 
exponent. 



Lemma 2.5. Let < p < 1. Then 



C2p 



< 



[r(i-p)] 2 „ 2 
r(2-2 P ) v 



Proof. By (^1]) and the definition of c ; 



A" 



|ntr p |(l-n)tr p dn = c 2 |tr 2p / u- p (l-u)- p du, 



u=0 



u=0 



and the result follows by evaluating the beta integral. 



□ 
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In particular, recalling T(l/2) = y^r, Lemmas 2.4 and 2.5 yield 

(2.7) m)\ < 



This proves (2.1) for p = 1, with ci < 47r, and thus by (|2.2| ) for every p < 1 with 
c p < (47r) p ; applying Lemma again, we obtain the finiteness of c p in (2.1) for all 
p < 2. Somewhat better numerical bounds are obtained for 1/2 < p < 1 by taking a 
geometric average between the cases p = 1/2 and p = 1: the inequality 

10(01 < (2t" 1/2 ) 2 ~ 2p (4^- 1 ) 2p - 1 = 2 2 V p - 1 t- p , t > 0, 

shows that c p < 2 2p 7r 2p ~ 1 , 1/2 < p < 1. In particular, we have c 3 / 4 < "V/87T, and thus, 
by Lemma |J, c 3/2 < M 1 ' 2 [r(l/4)] 2 < 186.4 < 187. 
Lemma 2.6. Le£ p > 1. T/ien 

c P+ i<2 p+1 C i+( 1 / P )p/( ;p _i). 

Proof. Assume that t > 2cJ v . Then, again using ((2] 

■ 1 



mm 



u=0 



(uty 



u=0 [(1 - u)t]P 



1 min 



[(1 - «)*]* 



, 1 I tin 



u=cl/*t-i [n(l-n)t 2 ]P 



1+(Vp) 



< 



< 



(1/2)p p 



p t p+i 
l+(i/ P ) t -(p+i) 



+ 2 



t 2p 



1/2 



i/i-t-i [tt(l - u)]p 



+ 



(i/2)p t 2p y u= £p t -i 



1/2 



u p du 



1 



eft 



^r 1 



-(p-i) 



= 2 p+1 c 1+(1/v) p r (p+1) . 
p p — 1 

1 lv 1 A? 

We have derived the desired bound for all t > 2c p . But also, for all < t < 2c p , we 
have 



2 p+i c i+(i/p) p t -(p+D > p 
p p — 1 ~ p — 1 



> 1 > 



so the estimate holds for all t > 0. 



□ 



Lemma ^1] completes the proof of finiteness of every c p in (2.1) (by induction), and 
of the estimate (2J3). The bound for C3/2 obtained above now shows (using Maple) that 
c 5/2 < 103215, which then gives c 7/2 < 197102280. 
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We can rewrite fl2.3|) as 

V(p+i) . _ 1 /„ / 1 



(p-i)(p+i) 



2 c y p exp( ' 



p K V2(p-l) 2(p+l) 
Hence, by induction, if p = n + | for a nonnegative integer n, then 



i/p < 2 n c^e (1/3)+(1/5) = C2 P , 



c 

where C := 2 _5//2 e 8//15 c^2 < 30.6 < 2 5 , using the above estimate of c 5 / 2 - Consequently, 
cj v < 2 P+5 when p = n + |. For general p > 3/2 we now use ( |2.2| ) with pi = p and 
p 2 = \p - |] + |, obtaining cj p < 2 P2+5 < 2 P+6 ; the case p < 3/2 follows from Q 
and the estimate < 33 < 2 6 . This completes the proof of fl2.4| ) and hence of 



Theorem 2.1. □ 



Remark 2.7. We used ( |1.1| ) in two different ways. In the first step we conditioned on 
the values of Y and Z, while in the inductive steps we conditioned on U. 

Remark 2.8. A variety of other bounds are possible. For example, if we begin with 



the inequality (P^Tj), use (|2.6|) , and proceed just as in the proof of Lemma |2.6| , we can 
easily derive the following result in the case t > 8ir: 

, \ . , . 32vr 2 / / t \ \ 32vr 2 lnt „ 

(2.8) \<f>(t)\ < _(]n(_)+ 2 ) < ^ for all t > 1.72. 

The result is trivial for 1.72 < t < 8ir, since then the bounds exceed unity 



Since Y has finite moments of all orders [recall (1.3)], the characteristic function 



is infinitely differentiable. Theorem 2.1 implies a rapid decrease of all derivatives, too. 
Theorem 2.9. For each real p > and integer k > 0, there is a constant c p ^ such that 

\(t> {k) (t)\ < c Pik \t\~ p for all t G R. 



Proof. The case k = is Theorem 2.1, and the case p = follows by \<)) {k) {t)\ < E\Y 



k 



The remaining cases follows from these cases by induction on k and the following calculus 
lemma. □ 

Lemma 2.10. Suppose that g is a complex-valued function on (0, oo) and that A,B,p > 
are such that \g(t)\ < At~ p and \g"(t)\ < B for all t > 0. Then \g'{t)\ < 2 v / ABt _p/2 . 
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Proof. Fix t > and let 9 = arg(</(i)). For s > t, 

Re(e- ie 5 '( S )) > Re( e -*V (t)) - |</( S ) - </(t)| > \g'(t)\ - B(s - t) 
and thus, integrating from t to t\ := t + (|(/(t)|/-B), 



Re(er l *( 5 (ti) - ff (t))) > / (b'WI - B(« - *)) ds 



Consequently, 



and the result follows. 



{t l -t)\g'{t)\-\B{t l -t) 2 = \g'{t)\ 2 /(2B). 



g'(t)\ 2 /(2B)<\g(t)\ + \g(t 1 )\<2At-P, 



□ 



In other words, the characteristic function <p belongs to the class S of infinitely dif- 
ferentiable functions that, together with all derivatives, decrease more rapidly than any 
power. (This is the important class of test functions for tempered distributions, intro- 



duced by Schwartz [14|; it is often called the class of rapidly decreasing C°° functions.) 



3 The limiting Quicksort density / and its derivatives 

We can now improve the result by Tan and Hadjicostas [15| on existence of a density / 
for Y. It is an immediate consequence of Theorem ^T, with p = and p = 2, say, 
that the characteristic function (j> is integrable over the real line. It is well-known — see, 
e.g., (3|, Theorem XV. 3. 3] — that this implies that Y has a bounded continuous density / 
given by the Fourier inversion formula 



(3.1) 



i 

2^ 



e~ ux 4>{t) dt, x G R. 



Moreover, using Theorem 2A with p = k + 2, we see that t k (j){t) is also integrable for 
each integer k > 0, which by a standard argument (cf. ||, Section XV. 4]) shows that / 
is infinitely smooth, with a kth. derivative (k > 0) given by 



(3.2) 



1 f°° 

f(k) {x)= / {-itfe-**<f>{t)db, 

2TT J t =-oo 

It follows further that the derivatives are bounded, with 



x G R. 



(3.3) 



sup 



|/( fc )(x)|<i- / \i\'- \oU)\<ll (/,•>()). 



t=— oo 



and these bounds in turn can be estimated using Theorem 2.1. Moreover, as is well 



known 14], [13, Theorem 7.4], an extension of this argument shows that the class S 
discussed at the end of Section is preserved by the Fourier transform, and thus The- 



orem 2.9 implies that / G S: 
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Theorem 3.1. The Quicksort limiting distribution has an infinitely differentiate den- 
sity function f . For each real p > and integer k > 0, there is a constant C p ^ such 
that 

\f {k) (x)\ < C Pik \x\- p for all x G R. □ 



For numerical bounds on /, we can use ( [3. 3D with k = and Theorem ^4] for several 
different p (in different intervals); for example, using p = 0, 1/2, 1, 3/2, 5/2, 7/2, and 
taking h = 4, t 2 = 4vr 2 , t 3 = (187/(4vr)) 2 , i 4 = 103215/187, t 5 = 197102280/103215, 

m<£- r m)\dt = - [°°m)\dt 

27T Jt^-oa IT J t=0 

1 f°° 

< - / min(l,2r 1/2 ,47rt- 1 ,187t~ 3/2 ,103215t" 5/2 ,197102280^ 7/2 )dt 
( 3 - 4 ) =l( f 1 dt + f 2 2t~ l l 2 dt+ 4nt- 1 dt+ t mr^dt 

7T \Jt=0 Jt=ti Jt=t 2 Jt=t 3 

+ / J 103215 1" 5/2 + / 197102280 £~ 7/2 

Jt=t 4 Vt=t 5 / 

< 18.2. 



Remark 3.2. We can do somewhat better by using the first bound in (|2.8| ) over the 



interval (103.18, 1984) instead of (as above) Theorem |27l| with p = 1, 3/2, 5/2, 7/2 over 
(103.18,i 3 ), (t 3 ,U), (t 4 ,t 5 ), (i 5 ,1984), respectively. This gives 

f(x) < 15.3. 

Similarly, (^~J) with k = 1 and the same estimates of 
1 f°° 1 



!/'(*)! < 



2vr 



i=— oo 



7T 



as in ( p.4[) yield 
dt < 3652.1, 



t=o 



which can be reduced to 2492.1 by proceeding as in Remark 3.2. The bound can be 
further improved to 2465.9 by using also p = 9/2. 

Somewhat better bounds are obtained by using more values of p in the estimates 
of the integrals, but the improvements obtained in this way seem to be slight. We 
summarize the bounds we have obtained. 

Theorem 3.3. The limiting Quicksort density function f satisfies max x fix) < 16 
and max x \f'{x)\ < 2466. □ 

The numerical bounds obtained here are far from sharp; examination of Figure 4 
of ]l5| suggests that max / < 1 and max \f'\ < 2. Our present technique cannot hope to 



produce a better bound on / than 4/ir > 1.27, since neither Lemma 2.3 nor (2.6) can 
improve on the bound \<fi(t)\ < 1 for \t\ < 4. Further, no technique based on (3.3) can 
hope to do better than the actual value of (2tt)~ 1 ff^_ ](/>(t)\ dt, which from cursory 
examination of Figure 6 of |15] appears to be about 2. 
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4 An integral equation for the density / 



Our estimates are readily used to justify rigorously the following functional equation. 

Theorem 4.1. The continuous limiting Quicksort density f satisfies (pointwise) the 
integral equation 

m=f I f{y)f( X - 9{u) - {l - U)y Y-dy d u, 

Ju=oJyeR, \ u J u 

where g(-) is as in (|1.2|). 

Proof. For each u with < u < 1, the random variable 
(4.1) uY + (1 - u)Z + g{u) 

[with notation as in ( |1.1[) 1 has the density function 

- 9{u) - (1 - u)z 



(4.2) 



fu{x) 



eR 



/(*)/ 



— dz, 

u 



where the integral converges for each x since, using Theorem 3.3, the integrand is 
bounded by / (z) (max f)/u < lQf(z)/u; dominated convergence using the continuity 
of / and the same bound shows further that f u is continuous. 

This argument yields the bound f u (x) < 16/ 'u, and since f u = f\- u by symmetry 
1]), we have f u {x) < 16/max(u, 1 



m 



u) < 32. This uniform bound, (Fl), and 
dominated convergence again imply that f f u (x) du is a continuous density for Y, and 
thus equals f(x) for every x. □ 

It was shown in jDJ that / is positive almost everywhere; we now can improve this 
by removing the qualifier "almost." 

Corollary 4.2. The continuous limiting Quicksort density function is everywhere pos- 
itive. 



Proof. We again use the notation (4.2) from the proof of Theorem 4.1. Fix x G R and 
u G (0, 1). Since / is almost everywhere positive flql , the integrand in (l4.2|) is positive 
almost everywhere. Therefore f u (x) > 0. Now we integrate over u G (0, 1) to conclude 
that f(x) > 0. □ 



Alternatively, Corollary |4.2| can be derived directly from Theorem |4.1| , without re- 
course to [p|. Indeed, if f(yo) > and uo G (0>1)> set x = yo + g(yo); then the 
integrand in the double integral for f(x) in Theorem [O] is postive for (u, y) equal to 
( M 0;?/o)) an d therefore, by continuity, also in some small neighborhood thereof. It fol- 
lows that /(yo + g(uo)) > 0. Since uq is arbitrary and the image of (0, 1) under g is 
(—(2 In 2 — 1), 1), an open interval containing the origin, Corollary [4.2| follows readily. 



Remark 4.3. In future work, we will use arguments similar to those of this paper, 
together with other arguments, to show that when one applies the method of successive 



substitutions to the integral equation in Theorem 4J, the iterates enjoy exponential- 
rate uniform convergence to /. This will settle an issue raised in the third paragraph of 



Section 3 in |15]. 
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